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1.  INTRODUCTION 


The  current  generation  of  interior  ballistic  codes  rely  on  turbulence  models  which 
incorporate  little  if  any  of  the  underlying  physics  of  the  gas— particle  motion  in  a  gun  tube. 
Indeed,  more  often  than  not,  the  fictitious  length  scales  and  adjustable  constants  employed 
are  based  on  experiments  which  have  no  relationship  to  the  flow  being  modeled.  The  result 
at  best  is  post  predictive  in  nature  and  no  claim  to  physical  representation  can  be  made. 
To  put  turbulence  sub— models  in  a  ballistic  setting  on  a  more  realistic  footing,  this  study 
was  undertaken.  A  new  mathematical  approach  to  the  problem  is  presented  and  results  for 
the  single  phase,  turbulent  case  are  given.  In  a  companion  paper  the  flow  in  a  two— phase 
jet  is  analyzed  and  compared  to  data  obtained  at  BRL. 

Turbulence  is  prevalent  in  practically  all  naturally  occurring  flow  processes  and  has 
challenged,  without  adequate  resolution,  the  best  scientific  minds  for  over  a  century.  The 
problem  is  especially  hard  due  to  the  very  large  number  of  coupled  degrees  of  freedom 
involved.  Hence  the  description  of  such  fluid  flows  cannot  proceed  via  usual  techniques  of 
approximating  the  system  by  a  linear  one  and  using  perturbation  techniques  around  the 
linearized  limit.  Despite  the  existence  of  an  enormous  body  of  theoretical  literature  on  the 
subject,  there  is  still  no  "standard"  theoretical  treatment  of  fully  developed  turbulence.  In 
our  previous  works  we  took  a  statistical  approach  to  this  problem  (Domokos, 
Domokos— Kovesi  and  Zoltani  1988),  following  the  pioneering  works  of  Martin,  Siggia,  Rose 
(1973)  and  De  Dominicis  and  PeliU  (1978).  In  this  approach  one  considers  the  classical 
equations  of  motion  of  a  fluid,  viz.  the  Navier-Stokes  equations,  perturbed  by  a  random 
force.  From  the  physical  point  of  view,  the  random  force  represents  the  fluctuations  in  the 
fluid;  its  presence  is  necessary  in  order  to  avoid  unstable,  non— turbulent  solutions  of  the 
equations  of  motion.  The  statistics  of  the  random  force,  in  turn,  generates  a  statistical 
distribution  of  the  components  of  the  flow  velocity.  The  correlation  functions  computed 
from  the  latter  distribution  can  be  compared  directly  with  the  results  of  measurements. 
An  application  of  the  technique  developed  by  these  authors  to  the  case  of  turbulent  channel 
flow  led  to  a  reasonable  agreement  with  the  available  experimental  data,  cf.  Burgett 
(1989).  The  approach  adopted  there  consisted  of  a  combination  of  analytical  and 
numerical  melhods:  an  analytic  approach  was  followed  and  numerical  computations  were 
used  only  for  the  computation  of  Fourier  transforms,  integrations,  etc.  While  this 
approach  led  to  a  considerable  reduction  of  computational  complexity,  it  became  clear  that 
a  straightforward  application  of  those  techniques  to  problems  like  fluid  flow  in  jets  would 
lead  to  an  unreasonably  large  amount  of  computing  time. 
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There  are  basically  two  features  of  the  approach  which  make  the  computations 
difficult.  First,  there  is  no  easy  way  available  to  satisfy  the  constraint  imposed  by  the 
equation  of  continuity  on  the  various  velocity  correlation  functions;  second,  the 
approximation  methods  used  are  essentially  perturbative  in  nature,  hence,  one  has  to  labor 
very  hard  in  order  to  achieve  a  reasonable  accuracy. 

In  this  work  we  describe  some  further  developments  in  the  statistical  theory  of  fully 
developed  turbulence  which  are  designed  to  avoid  the  difficulties  outlined  above.  First,  we 
describe  a  technique  by  means  of  which  the  equation  of  continuity  can  be  automatically 
satisfied.  (We  explicitly  describe  the  procedure  for  the  case  when  the  fluid  may  be 
regarded  an  incompressible  one,  although  the  method  can  be  generalized  to  arbitrary, 
compressible  flows.)  Second,  we  develop  a  variational  approach  to  the  computation  of 
correlation  functions.  This  enables  us  to  avoid  perturbative  approaches  altogether. 

The  paper  is  organized  as  follows.  In  the  next  section  we  briefly  review  the  statistical 
approach  to  turbulent  flows  as  described  in  Domokos,  Domokos— Kovesi  and  Zoltani  (198S). 
In  Sec.  3  we  introduce  vector  potentials  which  enable  us  to  satisfy  the  equation  of 
continuity  identically.  We  also  discuss  the  method  of  imposing  various  symmetry 
requirements  on  the  correlation  functions.  We  work  out  explicitly  the  constraints  imposed 
by  the  requirement  of  cylindrical  symmetry.  The  variational  approach  is  described  in  Sec. 
3.2.  A  sample  calculation  is  presented  in  Sec.  4:  some  correlation  functions  of  an 
axisymmetric  jet  are  calculated  by  means  of  the  techniqe  developed  here.  Sec.  5  contains  a 
discussion  of  the  results. 


2.  REVIEW  OF  THE  STATISTICAL  APPROACH 

Throughout  this  paper,  a  condensed  notation  is  used.  We  consider  a  class  of  systems 
described  by  a  classical  equation  of  motion  of  the  general  form: 

d  X  4-  F[X]  =  f(x).  (2.1) 

t 

Here  X  stands  for  an  element  of  the  vector  space  of  dynamical  variables.  F[X]  is  an 
autonomous  map  of  the  vector  space  upon  itself;  x  stands  for  both  spatial  coordinates  and 
the  time  variable,  x  =  (x,t).  Finally,  f(x)  represents  a  Gaussian  random  force  acting  upon 
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the  system;  it  plays  the  role  of  a  disordering  field,  driving  the  system  described  by  eq.  (2.1) 
away  from  a  non— chaotic  behavior.  Specifically,  if  the  system  is  an  incompressible  fluid 
described  by  the  Navier-Stokes  equations,  X  is  identified  with  the  velocity  field,  u.  The 
map  F[X]  reads: 


F[u]  =  (u  d  )u  +  d  p  —  yV2u  ,  (2.2) 

s 

with  the  velocity  field  satisfying  the  constraint,  V.u  =  0.  In  that  case,  the  quantity  p  (the 
pressure  divided  by  the  density),  is  not  an  independent  dynamical  variable:  it  can  be 
explicitly  expressed  in  terms  of  u.  Due  to  the  constraint  upon  u,  we  are  free  to  assume  that 
the  perturbing  random  force  is  solenoidal. 

The  quantity  playing  the  central  role  in  theories  of  this  type  is  the  generating 
functional  of  the  correlation  functions.  It  can  be  expressed  in  terms  of  the  probability 
distribution  of  the  random  force  as  follows.  Let  K  stand  for  the  correlation  operator  of  the 
random  force.  Then  the  generating  functional  of  the  correlation  functions  of  X  is  given  by 
the  functional  integral: 


Z[j]  =  |  DX  exp  -  [((ax  -  F),  K(a X  -  F))  +  (j,X)],  (2,3) 

where  (.,.)  stands  for  a  scalar  product  over  the  vector  space,  involving  integration  over 
space— time  variables  and  summation  over  tensor  indices,  whereas  j  is  an  arbitrary  function: 
the  functional  derivatives  of  Z  with  respect  to  j  give  the  correlation  functions.  The 
cumulants  are  generated  by  the  functional  W  =  — InZ.  Letting  f  to  be  a  white  noise,  viz. 
K(x[,x2)  =  k  £(ti  —  t2)<53(xi  —  X2)  (where  k  is  a  constant)  is  often  a  satisfactory  choice. 

The  functional  weight  of  integration,  DX,  is  proportional  to  the  (infinite) 
determinant, 


Det  ( d  -6F/5X).  (2.4) 

t 

It  was  shown  by  Domokos,  Domokos-Kovesy  and  Zoltani  (1988)  that  the  latter  can  be 
expressed  in  terms  of  a  functional  integral  over  Fadeev-Popov  ghosts  (Itzykson  and  Zuber 
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1980).  The  central  questions  are:  what  is  a  satisfactory  implementation  of  the  constraints 
imposed  upon  X  (in  the  case  of  an  incompressible  Quid  Qow,  the  constraint  V.u  =  0)  as  well 
as  the  development  of  suitable  approximation  techniques  for  the  computation  of  Z  (or  its 
functional  derivatives). 

3.  THE  MATHEMATICAL  APPROACH 

3.1  Incompressible  Flow:  Vector  Potentials  and  Symmetries.  An  incompressible  flow 
is  characterized  by  the  fact  that  the  velocity  field  is  solenoidal.  Any  solenoidal  vector  field 
is  obtained  as  the  curl  of  a  vector  potential.  Hence,  we  write,  u  =  V  *  A,  where  A  is  the 
vector  potential.  It  is  also  well  known  that,  given  the  velocity  field,  the  vector  potential  is 
determined  only  up  to  the  gradient  of  an  arbitrary  scalar  function:  A  and  A  +  V  H  give  rise 
to  the  the  same  velocity  field.  (This  is  called  the  gauge  freedom  in  electromagnetic  theory: 
we  adopt  the  same  terminology  here.)  The  scalar  function,  H  can  be  chosen  so  as  to 
simplify  the  problem  of  determining  the  vector  potential.  The  equation  satisfied  by  the 
vector  potential  is  obtained  by  substituting  the  relationship,  u  -  V  *  A  into  the 
Navier-Stokes  equation.  This  is  straight  forward  and  we  do  not  reproduce  the  result  here. 
Correspondingly,  in  the  statistical  theory  discussed  by  Domokos,  Domokos-Kovesi  and 
Zoltani  (1988)  and  briefiy  reviewed  in  the  preceding  section,  one  first  determines  the 
correlation  functions  of  the  vector  potential;  the  velocity  correlation  functions  are  then 
obtained  by  taking  the  curl  with  respect  to  every  argument. 

It  is  worth  remarking  that  this  way  of  satisfying  the  incompressibility  constraint  is 
not  the  only  pos  ability:  it  would  be  possible  to  take  the  constraint  into  account  by  means 
of  the  Fadeev— Popov  method,  see  Itzykson  and  Zuber  (1980).  However,  the  method 
described  here  is  simpler  and  leads  to  results  more  quickly  than  any  other  approach  we 
know  of.  This  is  due  to  the  gauge  freedom  just  described. 

In  what  follows,  we  use  a  gauge  such  that  the  component  of  the  vector  potential  along 
the  mean  Qow  vanishes.  ("Axial  gauge".)  On  denoting  the  mean  Qow  by  U,  the  condition 
to  be  satisfied  is:  U.A  =  0.  This  can  always  be  achieved.  In  fact,  let  B  be  a  vector 
potential  which  reproduces  the  velocity  field  ,  but  its  projection  onto  U  is  not  necessarily 
zero.  Then,  in  order  to  satisfy  U.A  =  0,  with  A  =  B  +  VH,  the  scalar  H  has  to  satisfy  the 
differential  equation, 

U.B  -f  U.VH  =  0.  (3.1.1) 
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There  remains  a  residual  gauge  freedom,  viz.  a  function  h  satisfying  U.Vh  =  0  can  always  be 
added  to  any  solution  of  (3.1). 


Any  symmetry  requirement  should  be  now  imposed  upon  the  correlation  functions  of 
the  vector  potentials  instead  of  the  correlation  functions  of  the  velocity  field  itself.  A 
symmetry  of  a  flow  means  that  the  correlation  functions  are  invariant  under  a  subgroup  of 
the  Euclidean  group.  For  instance,  homogeneous  turbulence  means  invariance  under 
translations,  hence,  the  n-point  correlation  function  of  the  vector  potential  depends  on  the 
n-i  coordinate  differences  only,  not  on  the  coordinates  themselves.  Likewise,  in  the  case  of 
isotropic  turbulence,  the  correlation  functions  of  order  2n  are  proportional  to  the  n— fold 
direct  product  of  the  metric  tensor  with  itself,  whereas  correlation  functions  of  odd  order 
vanish. 

Let  us  now  concentrate  on  the  two  point  correlation  function;  the  construction  of  the 
general  n-point  correlation  function  proceeds  in  a  similar  fashion.  For  the  sake  of 
definiteness,  we  choose  a  coordinate  system  such  that  its  third  axis  coincides  with  the 
direction  of  the  mean  flow,  U.  Due  to  the  fact  that  the  relationship  between  the  vector 
potential  and  the  velocity  field  is  a  linear  one,  a  Reynolds  decomposition  of  the  vector 
potential  leads  to  one  of  the  velocity  field.  Let  us  write, 

A  =  <A>  +  a, 

U  =  V  *  <A>, 

<a>  =  0, 

a’  =  V  *  a,  (3.1.2) 

where  u’  stands  for  the  fluctuating  part  of  the  velocity  field  and  <...>  denotes,  as  usual, 
the  expectation  value  of  a  quantity.  In  an  axial  gauge  we  have: 

A3  =  a3=0,  (3.1.3) 

so  that  only  the  components  of  the  vector  potential  transverse  to  the  mean  flow  are 
nonvanishing.  (These  components  are  denoted  by  subscripts/superscripts  in  capital  letters, 
A,  B,  etc.) 

The  two— point  correlation  function  of  the  vector  potential  a(x,t)  can  be  decomposed 
in  a  basis  of  second  rank  tensors  in  the  plane  perpendicular  to  the  mean  flow.  A  basis  of 
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such  tensors  consists  of  three  elements.  Correspondingly,  the  two— point  correlation 
function  can  be  written  as  follows: 


ZAB  =  <1A(ll'tl)  aB(l2'‘!)>  “  SAB  Zl  +  'AB  ^ 

+  Kx,axJb  +  x2ax,b)  Z„ 


(3.L4) 


where  6  and  t  are  the  Kronecker  and  Levi— Civita  tensors,  respectively.  The 

functions  Zj,  Z2  and  Z3  are  invariant  under  rotations  around  the  third  axis  and  they  are 
symmetric  functions  of  their  arguments,  (xptj)  and  (x2,t2)-  (The  antisymmetric  product  of 
the  two  coordinates  has  been  inserted  in  front  of  Z2  so  as  to  satisfy  the  permutation 
symmetry  of  the  correlation  function  with  all  three  invariant  functions  being  symmetric.) 

Thus  we  found  that  the  two  point  correlation  tensor  of  an  incompressible  fluid  has 
three  independent  components  only.  Cylindrical  symmetry  of  the  flow  further  reduces  the 
number  of  independent  elements  to  two.  In  fact,  cylindrical  symmetry  means  that  the 
correlation  function  is  invariant  under  rotations  around  the  third  axis;  the  tensors  e  and  t 
are  the  only  invariant  ones  under  such  rotations;  hence  Z3  =  0. 

We  now  give  the  explicit  expressions  of  the  velocity  correlation  functions  for  a  flow  of 
cylindrical  symmetry. 


G  = 
33 


G  -  <5AT3  di3<923  Z|, 

AB  AB 

G^2  =  —  $i3[xi^Z2  +  d2^Zi], 

RS  R  R 

^  +  X1  ^Ir.  +  X2  ^2r>]^2- 

Ho  HR 


(3.1.5) 


The  notation  used  in  this  and  subsequent  equations  is  the  following.  The  first 
subscript  denotes  the  argument  in  the  correlation  function;  the  second  subscript  or 
superscript  refers  to  the  vector  component.  The  summation  convention  is  used  throughout. 

Cylindrically  symmetric  turbulent  flows  of  incompressible  fluids  have  been,  of  course, 
discussed  previously,  cf.  Batchelor  (1946),  Chandrasekhar  (1950),  Trevino  (1982).  Where 
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they  overlap,  our  results  agree  with  those  of  these  works.  The  present  approach  is, 
however,  more  general,  since  it  can  be  applied  to  flows  of  arbitrary  symmetry  and  it  is 
more  explicit  than  the  treatment  of  Trevino  (1982). 

Finally,  we  just  mention  that  the  present  construction  can  be  generalized  in  a 
straight  forward  fashion  to  compressible  flows  where  both  the  velocity  field  and  the  density 
field  are  dynamical  variables.  In  that  case  one  has  to  introduce  two  vector  potentials  with 
a  correspondingly  enlarged  group  of  gauge  transformations;  this  problem  will  be  discussed 
elsewhere. 


3.2  Variational  Principles  for  the  Correlation  Functions.  In  this  Section  we  return 
to  the  condensed  notation  used  in  Sec.  2.  Consider  the  expression  of  the  generating 
functional  of  the  cumulants,  as  quoted  there.  We  have: 

w[j)  =  -  In  |  DX  exp  -  (((«  -  F),K(3  -  F))  +  (j,X)J,  (3.2.1) 

The  averages  of  the  various  quantities  are  given  by  functional  derivatives  of  W, 

M 

G(l)  =  <X(Xj)>  =  —  r  , 

*j(  x  l) 

£2W 

G ( 1 ,2)  =  <X(x,)  X(x2)>  =  -  ...  -  .  +  G(1)G(2)  ,  (3.2.2) 

<5j(xi)  <5  j(x2) 

etc.  Instead  of  the  functional  argument  j,  we  can  now  introduce  G(l)  as  a  functional 
argument,  by  means  of  a  Legendre  transformation  in  function  space.  Let  us  define  a  new 

functional, 


S,  =  W-(j,G,).  (3.2.3) 

One  readily  verifies  with  the  help  of  (3.2.1)  that  Si  is  a  functional  of  G(l),  its  first 
functional  derivative  being  given  by  the  expression: 
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(3.2.4) 


iS, 

io(T)  =  -j(I) 

The  functional  Si  can  be  determined  from  a  functional  differential  equation  obtained 
by  combining  eqs.  (3.2.1)  thru  (3.2.4).  This  gives: 

r  xo 

=  -  In  DX  exp  —  [((5  X  —  F),K(<9  X  —  Fl)  —  (-— ',X)].  (3.2.5) 

J  t  t  "  6  j 

The  differential  equation  can  be  solved  by  iteration;  one  usually  takes  as  a  zeroth 
approximation  the  functional: 

S,°  =  S[G ,]  =  ((t>  G,  -  F[G,]),K(fl  G,  -  F[G ,])  (3.2.6) 

X  X 

The  principal  advantage  of  the  functional  Sj  is  that  it  becomes  stationary  if  the 
external  ^source",  j,  is  put  equal  to  zero,  see  eq.  (3.2.4).  hence,  one  may  determine  Si  up 
to  a  certain  accuracy,  but  thereafter  a  calculation  of  Gj  does  not  have  to  rely  upon  any 
perturbative  approximation  scheme.  This  procedure  can  be  generalized  if  one  is  interested 
in  determining  the  higher  order  correlation  functions  as  well  as  the  average  of  the 
dynamical  variable  itself.  We  outline  the  procedure  for  a  scheme  aimed  at  determining  the 
functions  Gi  and  G2.  It  is  necessary  to  introduce  a  bilinear  source,  h(l,2),  in  addition  to  j, 
viz. 


W[j,h]  =  -  In  j  DX  exp  -  [((aX  -  F),K(0X  -  F))  +  (j,X)  +  (X,hX)] 

(3.2.7) 

Similarly  to  the  procedure  outlined  above,  one  now  performs  a  double  Legendre  transform 
in  order  to  obtain  a  functional,  Sa[G  1,62]: 

S2  =  W[j,h]-(Gi,j)-Tr(hG2),  (3.2.8) 
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where  the  trace  is  understood  in  the  operator  sense.  One  readily  verifies  the  relations, 


6  S2 
6G2 
SS_2_ 

SGt  __j' 

Thus,  S2  is  stationary  in  both  and  G2  in  the  limit  of  vanishing  sources.  The 
equation  obeyed  by  the  functional  S2  can  be  determined  in  the  same  way  as  eq.  (3.2.5).  We 
merely  quote  the  result: 


(3.2.9) 


S>-(G„§VTr(G,gi) 

=  -  In  |  DX  exp  -  [((3  X  -  F),K(3  X  -  F))  -  (^,X)  -  (X,^X)]. 

(3.2.10) 

Just  as  for  the  equation  obeyed  by  Si,  the  only  known  method  of  solving  (3.2.10)  is  by 
iteration.  We  give  here  the  first  approximation  to  S2: 

PS 

S21  =  S(G,1  +  Trl^-^-G,)  -  Tr(lnG,).  (3.2.11) 

(This  form  of  S21  is  obtained  by  changing  to  an  integration  variable,  Y  =  X  —  Gi  in  eq. 
(3.2.10);  thereafter  S[Y  +  G 1]  is  expanded  in  powers  of  Y  up  to  quadratic  terms.  The 
functional  differential  equation  is  then  solved  by  quadrature.) 

Functionals  of  the  type  (3.2.10)  were  first  used  by  De  Dominicis  and  Martin  (1964) 
and  Domokos  and  Suranyi  (1964);  the  differential  equations  obeyed  by  such  functionals  was 
obtained  by  Cornwall,  Jackiw  and  Tomboulis  (1974).  In  the  papers  just  quoted,  it  is  also 
shown  that  higher  order  correlation  functions  can  be  computed  by  the  application  of  either 
one  of  the  variational  principles  described  above.  In  this  work,  however,  we  concentrate 
upon  the  average  flow  (corresponding  to  Gj)  and  the  two  point  correlation  function.  For 
this  purpose,  the  use  of  the  variational  principle  based  on  the  functional  S2  is  best  suited. 
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4.  CALCULATION  OF  A  CYLINDRICALLY  SYMMETRIC  JET 


As  an  illustration  of  the  technique  developed  above,  we  now  perform  a  computation 
of  some  properties  of  a  free  jet,  with  a  cylindrically  symmetric  mean  flow.  We  use  the 
variational  principle  (3.2.8)  and  compute  the  functional  S2  in  the  first  approximation,  cf. 
eq.  (3.2.11).  The  power  of  the  variational  principle  lies  in  the  fact  that  one  can  use  a 
Rayleigh— Ritz  method  in  order  to  minimize  S2.  This  leads  to  a  considerable  reduction  in 
the  complexity  of  computation:  instead  of  solving  coupled  partial  differential  equations  for 
the  correlation  functions,  one  tries  to  guess  the  form  of  the  vector  potential  for  the  mean 
flow  and  of  the  functions  Zt  and  Z2  in  (3.1.5),  leaving  a  few  free  parameters.  On 
substituting  this  form  into  the  expression  of  S2,  the  problem  reduces  to  the  computation  of 
integrals  and  the  minimization  of  the  expression  so  obtained  as  a  function  of  the 
parameters.  The  success  of  such  an  approach  depends  on  finding  a  reasonable  funtional 
form  of  <A>  and  of  Zt  and  Z2.  This  can  be  accomplished  only  by  a  judicious  use  of  one’s 
physical  intuition  and  by  trial  and  error. 

In  this  calculation  we  want  to  simplify  the  calculation  as  much  as  possible:  we  guess 
some  simple  functional  forms  which  could  reproduce  the  main  features  of  the  experimental 
data. 


Let  us  examine  now  the  Navier— Stokes  operator  entering  the  expression  of  the 
functional  S2: 


N  [u]  =  d  u  -  vVH  +  (ud  )u*  +  dp  (4.1) 

t  s 

We  are  interested  in  fully  developed  turbulence:  that  means  that  it  is  described  by  a 
stationary  ensemble.  Hence  the  term  containing  the  time  derivative  can  be  omitted  from 
(4.1).  Next  we  notice  that  if  one  introduces  dimensionless  quantities  (as  we  shall  do  in 
what  follows),  the  coefficient  of  the  viscous  term,  V2u,  becomes  proportional  to  1/Re,  where 
Re  stands  for  the  Reynolds  number.  Hence,  for  flows  of  large  Reynolds  numbers  (Re  a  104, 
say),  it  is  safe  to  omit  the  viscous  term  too,  as  long  as  we  are  interested  in  the  behavior  of 
the  fluid  on  scales  substantially  larger  than  the  dissipation  scale.  Thus,  we  are  going  to 
work  with  the  truncated  Navier— Stokes  operator, 
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with 


(4.2) 


n^u]  +  dp, 
n*[u]  =  u  d  u*. 


The  next  task  is  to  eliminate  the  pressure  from  the  generating  functional.  To  this 
end,  we  regard  the  pressure  as  one  of  the  fluctuating  dynamical  variables  and  carry  out  the 
functional  integration  over  it  explicitly.  Using  an  operator  notation  as  before,  we  have  to 
compute  a  functional  integral  of  the  form: 


r 

j 


Dp  exp— j[n  +  Vp]  K  [n  +  Vpj. 


(4.3) 


This  is  a  standard  Gaussian  integral.  The  result  of  the  computation,  to  be  inserted  into 
the  expression  of  the  functional  S2,  is  the  following: 


with 


f 

J 


T 

Dp  exp-j[n  +  Vp]  K  [n  +  Vp]  =  exp— jn  K  n, 


=  K  -  (VK)'  [VKV  +  a2] 


VK, 


(4.4) 


(4.5) 


where  '  denotes  the  transpose  of  an  operator. 


In  the  last  equation  we  inserted  a  constant  a2  in  order  to  make  the  inverse  of  the 
operator  VKV  well  defined  on  large  length  scales  (equivalently,  at  small  wave  numbers). 
T  T 

Evidently,  K  is  the  transversepart  of  K,  VK  =  0.  We  now  have  to  make  a  physical 
assumption  about  the  correlation  operator,  K.  We  argue  that  the  stirring  force  should 
point  in  the  direction  of  the  mean  flow:  this  is  both  intuitively  plausible  and  it  is  the 
simplest  wray  of  satisfying  the  requirement  of  cylindrical  symmetry  of  the  problem. 
Otherwise,  we  assume  the  the  correlation  to  be  of  short  range.  Hence  we  take,  (cf.  Sec.  2): 
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K  (x,x’)  =  k  6  S  63(x  —  x’). 
ij  '  i3  j3  V  ' 


(4.6) 


In  order  to  obtain  a  finite  result  even  in  the  limit  a  ->  0,  we  have  to  assume  that  k  is 
proportional  to  1  /a.  With  this,  the  transverse  part  of  K  becomes: 


K 


r  .  ,  —a |z  —  z’ I  A  A. 

i,i’  =ae  £2(x  —  x’  )5  S  , 

ijV  i3  j3 


(4.7) 


where  a  is  a  constant  independent  of  a.  At  the  end  of  the  calculation  one  may  take  a  ->  0, 
although  there  is  some  evidence  that  a  finite  a  (just  as  in  an  Ornstein— Zernike  process) 
gives  slightly  better  results,  cf.  Burgett  (1989),  Ch.6. 

We  are  now  ready  to  compute  the  functional  S21,  eq.  (3.2.11).  Clearly,  the  first  term 

reads: 


S[U]  =  /d3xd3x;  n  [U(x)]K^.  .(x,x’)n^[U(x’)] 


(4.8) 


On  expressing  the  average  velocity  as  the  curl  of  its  vector  potential,  U  =  V  *  A  and 
using  (4.6),  we  obtain  in  the  axial  gauge: 


q  t A2  a  a  f  AS  RM 

S  =  a  /  d2x  dzi  dz2  e  [  —  e  o  A  t  do  A 

1  3  S  ARM 

R  S 

+  R(£  aRA^  I-]  > 

3  R  a  Z\  Z2 


(4.9) 


where  the  second  factor  in  square  brackets  has  the  same  structure  as  the  first  one,  but  it  is 
evaluated  at  Z2.  We  now  notice  that  the  second  term  in  both  square  brackets  in  the  last 
equation  is  a  total  derivative.  Consequently,  upon  integration  by  parts,  its  contribution 
becomes  proportional  to  a  and  hence  it  is  small  for  small  values  of  that  parameter. 
Therefore,  that  term  can  be  omitted  without  substantially  affecting  the  iesults  and  we  shall 
do  so  in  this  paper. 
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We  next  compute  the  functional  derivatives  of  (4.9)  in  order  to  generate  the  second 
term  of  S21  in  eq.  (3.2.11).  This  is  a  straight  forward  procedure  and  we  merely  quote  the 
result. 


<52S  -a  Z[  -  z2  NA  RB 

7—7 - r~7T  ; - r=2ae  {2ae(z1-z2)£  c  d  d  6\x i-x2) 

5A  (x^zi)  <5A  (x2,z2)  R  N 


B' 

NA  RN1  PB  QSC  Q  4 

+  a!‘  e  VrV  ‘  VqV(x,_X!)>' 


(4.10). 


In  this  equation  and  from  now  on,  xi,  etc.  stand  for  the  transverse  components  of  the 
position  vector,  whereas  zt,  etc.  denote  the  longitudinal  component  of  the  same  vector. 
The  symbol  e(z)  stands  for  the  sign  function,  e(z)  =  z/|z|.  Unlike  in  the  expression  of  S, 
we  now  cannot  omit  terms  of  0(a)  or  0(a2),  for  there  are  no  terms  of  0(1)  present  in  eq. 
(4.10).  However,  we  notice  that  the  quantity  a  can  be  scaled  out  of  the  expression  S21.  In 
fact,  on  looking  at  the  structure  of  eqs.  (3.1.4),  (3.2.11),  (4.9)  and  (4. 10), we  realize  that 

upon  the  rescaling,  z  -*  (l/a)z,  A  -*  (l/o^)A,  and  Z  -*  (l/a3)Z  ,  the  expression  of  S21 

AB  AB 

is  multiplied  by  an  overall  factor  1/a4.  This,  however,  can  be  absorbed  into  the  constant 
entering  the  expression  of  the  force  correlation  function;  at  any  rate,  the  location  of  the 
stationary  point  of  the  functional  S2l  is  independent  of  the  value  of  a.  (We  remark, 
however,  that  higher  order  approximations  to  S2  do  not  have  this  scaling  property.) 

The  reader  certainly  notices  that  our  approach  has  been  quite  general  and  the  results 
obtained  so  far  are  applicable  to  almost  any  flow  geometry.  We  now  specialize  to  the  case 
of  a  cylindrically  symmetric  free  jet.  Such  jets  have  a  characteristic  velocity,  namely  the 
mean  flow  velocity  at  the  centerline  of  the  jet  exit,  Uom,  and  a  characteristic  length,  d, 
namely  the  diameter  of  the  jet  at  the  exit.  From  now  on,  as  is  customary,  we  measure  all 
velocities  and  distances  in  these  units.  However,  in  order  to  keep  the  notation  simple,  this 
is  not  indicated  explicitly  in  the  formulae.  (For  instance,  a  velocity  component  denoted  by 
U  is  understood  to  mean  U/Uom,  or  an  axial  distance,  z,  is  to  be  interpreted  as  z/d.) 

A  cylindrically  symmetric  mean  flow  can  be  described  in  the  axial  gauge  by  means  of 
a  vector  potential  which  has  a  tangential  component  only.  Specifically,  in  cylindrical 
coordinates  (xl  =  r  cos<p,  x2  =  r  sinip),  we  choose  a  trial  function  for  the  vector  potential  of 
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the  form: 


A  =  *U(0,z)  R(z)2  [1  -exp(-r2/R(z)2)],  (4.11) 

with  the  other  two  components  vanishing.  Here  U(0,z)  stands  for  the  mean  velocity  along 
the  jet  axis;  the  quantity  R(z)  characterizes  the  radius  of  the  jet  at  distance  z  from  the 
exit.  The  choice  of  such  a  functional  form  is  motivated  by  its  simplicity  and  by  various  fits 
and  approximate  calculations  of  free  jets;  see,  e.g.  the  classic  text  of  Hinze  (1975).  This 
vector  potential  gives  rise  to  an  axial  and  a  radial  component  of  the  mean  velocity  of  the 
form: 


Uz(r,z)  =  U(0,z)  e 


— r2/R(z)2 


„  ,  ,  dlnU(O.z)  -r2/R(z)\ 

Ur(r,z)  =  { — — iR(z)!  (1  -  e  ) 

,  ,dR,  ,  r2  — r2/R(z)2 
+  R(z)-[1-(1  +  -)e  U]}. 


(4.12) 


Next  we  impose  the  requirement  of  cylindrical  symmetry  on  the  two-point 

correlation  function,  Z  ,  cf.  eq.  (3.1.4): 

AB  v  ' 


(4.13) 


Zi  + 


AB 

e 


z2. 


In  order  to  simplify  the  computation,  we  arbitrarily  set  Z2  =  0.  There  is  no  strong  physical 
motivation  for  this  choice;  it  simplifies  the  calculations  to  a  considerable  extent,  albeit  at 
the  cost  of  some  loss  of  accuracy.  We  assume  the  following  functional  form  for  Zp 


=  A  e 


..  A  A  .  A  A 

-f(Z! -Z2)2 -g(Xi  -X2  )(xi  X2  ) 


«  (U(0,zt)U(0,z2)jd  (1  +  BM)  e 


(4.14) 


In  this  equation,  the  parameters  A,  B,  d,  6,  f,  and  g  are  treated  as  variational  parameters, 
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although  later  we  found  that  the  results  are  rather  insensitive  to  the  choice  of  d  and  we 
ended  up  using  a  plausible  value,  d  =  1/2.  The  expression  M  is  defined  as  follows: 


M  - 


i( 


A  A 
x  1  x  i 

R  ( z  i) 2 


+ 


A  A 

*2  X2  ' 

R(  z2) 2 


(4.15) 


The  radius  of  the  jet  was  assumed  to  be  a  linear  function  of  z,  R(z)  =  a  +  bz,  with  the 
parameters  a  and  b  being  also  varied.  We  did  not,  however,  vary  the  functional  form  of  the 
mean  velocity  on  the  jet  axis:  although  this  is  possible  in  principle,  we  used  a  fixed  form  in 
order  to  simplify  the  calculational  task.  We  used  the  following  simple  modification  of 
Spalding’s  formula  (Hinze  1975),  which  is  quite  accurate  for  z  >  2: 


U(0,z) 


1.35 


1  +  0.038z 


(4.16) 


The  next  task  is  to  insert  the  expressions  of  the  vector  potential  of  the  average  flow 
and  of  the  two  point  correlation  function  into  the  functional  S21  and  to  minimize  that 
expression  with  respect  to  the  parameters.  The  operations  are  elementary,  but  extremely 
tedious  to  carry  out  by  hand,  although  this  is  not  impossible.  They  are,  however,  easily 
carried  out  by  a  symbolic  manipulation  program.  We  carried  out  the  differentiations  and  a 
part  of  the  integrations  by  using  the  program  "Maple",  (Char  et  al.  1985).  The  search  for 
an  extremum  of  the  functional  was  done  numerically.  All  computations  were  done  on  a 
Sun  3/160  computer.  In  order  to  carry  out  the  search  for  the  extremum,  one  needs  starting 
values  of  the  parameters.  These  were  obtained  by  comparing  Uz(r,z),  eq.  (4.12),  to  a  few 
data  points  from  Zoltani  and  Bicen  (1990a);  we  also  chose,  arbitrarily,  f  =  g  =  6  =  1  as 
starting  values.  (The  values  of  f  and  g  cannot  be  read  off  directly  from  the  data.) 

The  search  resulted  in  the  following  values  of  the  parameters: 

a  =  0.25,  b  =  0.076; 

A  =  0.0013,  B  =  2.3 
6 =  1.79, 

f#5.6,  g»3.7  (4.17) 
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In  the  approximation  used,  the  value  of  the  functional  S2  changes  rather  slowly  vith 
variations  of  the  parameters  f  and  g.  This  is  due  to  the  fact  that  in  this  approximation,  the 
interaction  between  the  fluctuations  and  the  mean  flow  is  taken  into  account,  but  the 
mutual  interaction  of  the  fluctuations  is  neglected.  (The  functional  S21  is  almost  local:  it 
contains  6  functions  and  derivatives  of  6  functions  of  finite  order.)  Hence,  the  values  of  f 
and  g  cannot  be  determined  efficiently.  This  is  not  a  serious  deficiency,  however:  most 
measurements  determine  correlation  functions  at  coincident  arguments,  which  are  rather 
insensitive  to  these  parameters. 

The  correlation  functions  of  velocity  fluctuations  are  determined  by  inserting  the 
expression  of  Z  into  eq.  (3.1.5).  We  summarize  the  results  by  listing  the  expressions  of  the 
axial  component  of  the  mean  flow  and  of  the  correlation  functions  at  coincident  arguments. 

G„  =  G22  =  2fA  U2(0,z)  e-*2{l  +  BX 2  +  (l/2f)[d(z)2 
+  X2  B  d(z)  (d(z)  -  (B— 1)/B  b/R(z)) 

+  X4  B  b/R  (d(z)  +  (1  -  2B)/4B  b/R(z)) 

+  X5B/ 4  (b/R(z))2]}, 

G33  =  4gA  Uz(0,z)  e“X2  {1  +  x2(B  -  (2B  -  l)/4R(z)2g] 

+  X4B/4R(z)2g}, 

G13  =  G23  =  A  X  e_*2  U(0,z)  {  (1  -  B)  d(z) 

+  x2  [B  d(z)  +  (1  -  2B)  b/R]  4-  x4  bB/R(z)}, 

Uz(r,z)  =  U(0,z)  e_X2.  (4.18) 

Here,  x2  =  Ar2/R(z)2and  d(z)  =  d/dz(U(0,z))^. 

The  results  listed  in  eq.  (4.18)  are  evaluated  using  the  parameter  values  (4.17)  and 
compared  with  data  taken  from  Zoltani  and  Bicen  (1990a),  (1990b)  in  Figures  1  thru  4. 
(The  experimental  circumstances  there  were  not  exactly  identical.  In  particular,  the 
Reynolds  numbers  differ  by  approximately  30%  at  the  jet  exit.  However,  at  large  Reynolds 
numbers  the  profile  of  the  mean  flow  and  the  correlation  functions  should  be  insensitive  to 
the  precise  value  of  Re,  cf.  the  discussion  at  the  beginning  of  this  Section.  The  data  bear 
out  this  conclusion.)  In  the  figures  we  use  conventional  notation.  The  correspondence 
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between  the  notation  used  in  this  paper  (better  suited  for  theoretical  calculations)  and  the 
conventional  one,  as  used  in  Zoltani  and  Bicen  (1990a),  (1990b)  is  the  following.  (All 
arguments  are  suppressed,  so,  for  instance,  Gn  stands  for  Gn(x,x).  This  does  not  lead  to 
any  ambiguity.) 


G33  = 
G 11  = 

G22  = 
G 13  = 


<u2> 

U  0  ID  2 ' 

<v2> 
Uom2  ’ 

<W2> 

U  0  ID  2  ’ 

<uv> 
Uom2  ’ 


u,  =  u 


(4.19) 


5.  DISCUSSION 

A  look  at  Figures  1  thru  4  shows  that  overall,  the  simple  trial  functions  used  in  our 
variational  calculation,  with  the  parameters  determined  by  extremizing  S2,  are  in  a 
reasonable  agreement  with  the  data.  (The  agreement  with  the  data  is  better  at  large 
values  of  z.  This  is  understandable  on  physical  grounds:  at  larger  values  of  z,  the  statistics 
of  the  flow  is  closer  to  the  stationary  ensemble  we  have  been  working  with.)  In  particular, 
our  assumption  of  cylindrically  symmetric  correlation  functions  appears  to  be  well  justified. 
This  is  expected  on  physical  grounds:  to  a  very  good  approximation,  the  mean  flow  is 
cylindrically  symmetric.  Although  this  does  not  necessary  imply  that  the  correlation 
functions  themselves  should  possess  the  symmetry,  one  expects  violations  of  this  symmetry 
to  occur  on  relatively  short  time  scales.  Likewise,  the  assumption  about  "transverse 
scaling",  namely  that  the  radial  dependence  of  both  the  mean  flow  and  the  correlation 
functions  occurs  only  in  the  scale  invariant  combination,  x 2  =  <fr2/R(z)2  appears  to  be 
reasonably  well  satisfied. 

A  notable  exception  is  the  correlation  function  G 13,  which  is  predicted  to  be  much  too 
small  compared  to  the  data.  Various  modifications  of  the  form  of  the  trial  functions  have 
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Figure  1.  Predicted  and  Measured  Mean  Velocity  Profiles  at  fz/d)  =  5.  The  Values  Were 
Normalized  With  the  Centerline  Jet  Exit  Velocity. 
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Figure  3.  Predicted  and  Measured  Autocorrelation  Function.  <v  >  in  the  Radial 
Direction  at  (z/'d)  =  5.  The  Values  Were  Normalized  With  the  Square  of  the  Centerline  Jet 
Exit  Velocity. 
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1  Allegheny  Ballistics  Laboratory 
Propulsion  Technology  Department 
Hercules  Aerospace  Company 
ATTN:  Mr.  Thomas  F.  Farabaugh 
P.O.  Box  210 
Rocket  Center,  WV  26726 

1  MBR  Research  Inc. 
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Cdr,  CSTA 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes. 
Your  comments/answers  to  the  items/questions  below  will  aid  us  in  our  efforts. 

1 .  BRL  Report  Number _ BRL-TR-3192 _ Date  of  Report  JANUARY  1991 

2.  Date  Report  Received _ 

3.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest 

for  which  the  report  will  be  used.) _ 


4.  Specifically,  how  is  the  report  being  used?  (Information  source,  design  da»a,  procedure,  source 
of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars 
saved,  operating  costs  avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate.  _ 


Yi.  Gerferal  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate 
changes*  W  organization,  technical  content,  format,  etc.) _ 


-  Y 

‘ 

Name 


CURRENT  Organization 

ADDRESS _ 

Address 


City,  f  (ate,  Zip  Code 

7  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  New  or  Correct 
■'d'Jress  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


Name 


OLD  Organization 

ADDRESS  _ 

Address 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  mail.) 
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Director 

U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  SLCBR  DD-T 

Aberdeen  Proving  Ground,  MD  21(>i  -50fS6 
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BUSINESS  REPLY  MAIL 

FIRST  CLASS  PERMIT  No  0001,  APG,  MD 


POSTAGE  WILL  BE  PAID  BY  ADDRESSEE 


Director 

U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  SLCBR-DD-T 

Aberdeen  Proving  Ground,  MD  21005-9989 


FOLD  HERE 


